Effects of diversity on thermal niche variation in bird communities under climate change

Climate change alters ecological communities by affecting individual species and interactions between species. However, the impacts of climate change may be buffered by community diversity: diverse communities may be more resistant to climate-driven perturbations than simple communities. Here, we assess how diversity influences long-term thermal niche variation in communities under climate change. We use 50-year continental-scale data on bird communities during breeding and non-breeding seasons to quantify the communities’ thermal variability. Thermal variability is measured as the temporal change in the community’s average thermal niche and it indicates community’s response to climate change. Then, we study how the thermal variability varies as a function of taxonomic, functional, and evolutionary diversity using linear models. We find that communities with low thermal niche variation have higher functional diversity, with this pattern being measurable in the non-breeding but not in the breeding season. Given the expected increase in seasonal variation in the future climate, the differences in bird communities’ thermal variability between breeding and non-breeding seasons may grow wider. Importantly, our results suggest that functionally diverse wildlife communities can mitigate effects of climate change by hindering changes in thermal niche variability, which underscores the importance of addressing the climate and biodiversity crises together.

. Variation in community diversity among survey routes/sites within grid cells (N = 65). Panels on the left represent taxonomic, functional, and evolutionary diversity measures in the breeding season. Panels on the right represent taxonomic, functional, and evolutionary diversity measures in the non-breeding season. Here, taxonomic diversity is represented as the number of species, functional diversity as functional richness, and evolutionary diversity as phylogenetic diversity (PD). For visualization, outlier values are omitted. Figure S2. Variation in the mean phylogenetic diversity (PD) values in each grid cell calculated based on 25 replicated phylogenetic trees in breeding and non-breeding seasons. Phylogenetic diversity is measured as the sum of phylogenetic tree branch lengths in the phylogenetic tree of the community and then averaged over the survey routes/sites within the grid cell. Each box plot represents one grid cell (N = 65). Figure S3. Pairwise correlations of taxonomic, functional, and evolutionary diversity measures quantified with Pearson's correlation coefficients. The upper and lower panels describe the diversity measure correlations within grid cells for breeding and non-breeding season data, respectively. Diet specialization = Community weighted mean of diet diversity, Habitat specialization = Community weighted mean of habitat niche breadth, Temperature specialization = Community weighted mean of thermal niche breadth, Vertebrates in diet = Community weighted mean of vertebrates in diet. For description of diversity measures, see 'Methods' in the main text.  (1966)(1967)(1968)(1969)(1970)(1971)(1972)(1973)(1974)(1975)(1976) and in the full dataset  in breeding and non-breeding seasons.  Figure S4. Temporal trends in climate across the study area between 1966 and 2016. As a proxy for measuring regional climate change, we used the observed monthly temperature anomalies in 5° × 5° grid cell data from Earth System Research Laboratory (Jones et al., 2012). Panels a) and b) show the linear temperature trend in the breeding and non-breeding season data, respectively. In both panels, hexagon color gradient indicates the number of observations (here, sampling routes/sites in a year) and the dashed line indicates the linear fit of the regression model: temperature anomaly ~ year. Figure S5. Temporal variation in diversity measures in the breeding season bird communities. Each panel shows the temporal variation in one diversity measure (indicated on the y-axis). In all panels, hexagon color gradient indicates the number of observations (here, sampling routes in a year) and the dashed line indicates the linear fit of the regression model: diversity measure ~ year. In panels d-g, CWM stands for communityweighted mean of the particular trait.   Red circles illustrate the grid cells that had data of the non-breeding season bird communities, but not of the breeding season bird communities. On the contrary, blue circles illustrate the grid cells that had data of the breeding season bird communities, but not of the non-breeding season bird communities.

Supplementary results
Among all the diversity-thermal variability models (response variable: |∆CTI|), models including community diversity measures performed better than the null model in the non-breeding but not in the breeding season (Table S3). In the breeding season, the null model performed best, while the second best model included the functional dispersion. In the non-breeding season, the diversitythermal variability models that performed the best included the community weighted mean of vertebrates in diet, while the second best model included the community weighted mean of habitat niche breadth. In the non-breeding season, |∆CTI| was smaller when diversity was higher (Table  S4). We mainly did not find differences in the diversity-thermal variability relationships between breeding and non-breeding seasons (Table S5). More specifically, the interaction term was not significant in any of the models, except in the case of functional richness. The relationship between thermal variability and functional richness was stronger in the non-breeding season.
We inspected visually the best models for each season and found the model residuals normal and variances homogenetic in case of the best model fitted to non-breeding season data (Figures S4-S5). Table S3. Output of linear models on the diversity-thermal variability relationships (response variable: |∆CTI|) between 1966 and 2016. Each row represents one hypothesis and one linear model with the same response variable and different explanatory variables. All models include a Gaussian spatial correlation structure and are weighted by the grid cell-specific standard errors obtained from the models CTI ~ Year. The best model(s) were identified using Akaike Information Criteria for small samples (AICc). The rows for each season are ordered from best to worst model according to ∆AICc, starting with the best supported model (∆AICc=0). We identified as the best model(s) those models for which ∆AICc ≤ 2. Pseudo R 2 for gls-models are represented to indicate models' explanatory power. The coefficient and standard error of diversity measure variables in each model are shown in separate columns. Similarly, the coefficient and p-value of temperature change variable in each model are shown in separate columns.  Table S4. Linear model outputs for the best diversity-thermal variability relationships (response variable: |∆CTI|) in the breeding (B) and non-breeding (NB) seasons (for all models, see Table S3). All models include a Gaussian spatial correlation structure and are weighted by the grid cellspecific standard errors obtained from the models CTI ~ Year.   Figure S9. Normality of residuals in the best models (response variable: |∆CTI|) in the breeding (B) and non-breeding seasons (NB) (for the ranking of all models, see Table S3). Each panel represents one of the models listed in Table S4. Figure S10. Variance homogeneity of residuals of the best models (response variable: |∆CTI|) in the breeding (B) and non-breeding (NB) seasons (for the ranking of all models, see Table S3). Each panel represents one of the models listed in Table S4.

Sensitivity analyses
To account for the potential variation in temporal trend strengths of CTI and in the development of data collection methods, we conducted a sensitivity analysis of thermal variability models using subsets of the full datasets. By 1980, use of optical tools for bird observations became a standard, potentially increasing the detection probabilities. We refitted the diversity-thermal variability models using thermal variability measure and community data between years 1980 and 2016. We found that |∆CTI| did not differ between the full and subset data analyses, such that the means across the grid cells in the breeding season were 0.009 and 0.009 and in the non-breeding season were 0.036 and 0.036, in the full and subset data analyses respectively. Qualitatively the general diversity-thermal variability results remained parallel in the sensitivity analysis as in the full data analysis (compare Figures S11 and Figure 4 in the main text). However, the ranking of the models differed, such that in the breeding season, the temperature model was ranked the best and in the non-breeding season, the evolutionary and taxonomic diversity models were ranked above other models (Table S9). However, the functional richness model in the non-breeding season had the steepest diversity-thermal variability slope but was not top-ranked due to high variation (Table S9).
We also tested for the relationship between temperature change and thermal variability in each season and found no significant relationships ( Figure S12). Table S9. Output of linear models on the diversity-thermal variability relationships (response variable: |∆CTI|) for data between years 1980 and 2016. Each row represents one hypothesis and one linear model with the same response variable and different explanatory variables. All models include a Gaussian spatial correlation structure and are weighted by the grid cell-specific standard errors obtained from the models CTI ~ Year. The best model(s) were identified using Akaike Information Criteria for small samples. The rows for each season are ordered from best to worst model according to ∆AICc, starting with the best supported model (∆AICc=0). We identified as the best model(s) those models for which ∆AICc ≤ 2. Pseudo R 2 for gls-models are represented to indicate models' explanatory power. The coefficient and standard error of diversity measure variables in each model are shown in separate columns.  Figure S11. Examples of diversity-thermal variability relationships (response variable: |∆CTI|) in a)-c) breeding and d)-f) non-breeding seasons for data between years 1980 and 2016. The linear relationships with metrics of taxonomic, functional and evolutionary diversity for each season are illustrated. Thereby, panels a) and d) illustrate the temporal change in community temperature index (|∆CTI|) as a function of the number of species, panels b) and e) illustrate |∆CTI| as a function of the functional dispersion, and panels c) and f) illustrate |∆CTI| as a function of the phylogenetic diversity within the community. The model explanatory power and statistical significance are labelled at the top of each panel. Note that the effect of temperature trend is not accounted for in the model structure of the illustrated relationships (see Table S7).